home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / malign.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  9.3 KB  |  380 lines  |  [TEXT/R*ch]
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <ctype.h>
  4. #include <stdlib.h>
  5. #include "clustalw.h"
  6.  
  7. /*
  8.  *       Prototypes
  9.  */
  10. extern  void    *ckalloc(size_t);
  11. extern  void    ckfree(void *);
  12. extern  int      prfalign(int *, int *);
  13. extern  void    calc_seq_weights(int, int);
  14. extern  int    calc_similarities(int);
  15. extern  void    create_sets(int, int);
  16. extern  void    aln_score(void);
  17. extern  void    clear_tree(treeptr);
  18. extern  int    read_tree(char *, int, int);
  19.  
  20. int malign(int istart);
  21. int palign1(void);
  22. int palign2(void);
  23. double countid(int s1, int s2);
  24.  
  25. /*
  26.  *       Global Variables
  27.  */
  28.  
  29. extern double  **tmat;
  30. extern int     debug;
  31. extern int     max_aa;
  32. extern int     nseqs;
  33. extern int     profile1_nseqs;
  34. extern int     nsets;
  35. extern int     **sets;
  36. extern int     divergence_cutoff;
  37. extern int     *seq_weight;
  38. extern int     output_order, *output_index;
  39. extern Boolean distance_tree;
  40. extern char    phylip_phy_tree_name[];
  41. extern char    seqname[],treename[];
  42. extern int     *seqlen_array;
  43. extern char    **seq_array;
  44.  
  45. int malign(int istart) /* full progressive alignment*/
  46. {
  47.    static     int *aligned;
  48.    static     int *group;
  49.    static     int ix;
  50.  
  51.    double     *maxid, max;
  52.    int         val,i,j,k,ns,len,set,status,iseq;
  53.    int         pos,entries,ptr;
  54.    int        score = 0;
  55.  
  56.  
  57.    fprintf(stdout,"\nStart of Multiple Alignment\n");
  58.  
  59.    seq_weight = (int *) ckalloc( (nseqs) * sizeof(int) );
  60.  
  61. /* get the phylogenetic tree from *.ph */
  62.  
  63.    if (nseqs > 3) 
  64.      {
  65.        status = read_tree(phylip_phy_tree_name, 0, nseqs);
  66.        if (status == 0) return(0);
  67.      }
  68.  
  69. /* calculate sequence weights according to branch lengths of the tree -
  70.    weights in global variable seq_weight normalised to sum to 100 */
  71.  
  72.    calc_seq_weights(0, nseqs);
  73.  
  74. /* recalculate tmat matrix as percent similarity matrix */
  75.  
  76.    status = calc_similarities(nseqs);
  77.    if (status == 0) return(0);
  78.  
  79. /* for each sequence, find the most closely related sequence */
  80.  
  81.    maxid = (double *)ckalloc( (nseqs+1) * sizeof (double));
  82.    for (i=1;i<=nseqs;i++)
  83.      {
  84.          maxid[i] = 0.0;
  85.          for (j=1;j<=nseqs;j++) 
  86.            if (maxid[i] < tmat[i][j]) maxid[i] = tmat[i][j];
  87.      }
  88.  
  89. /* group the sequences according to their relative divergence */
  90.  
  91.    if (istart == 0)
  92.      {
  93.         sets = (int **) ckalloc( (nseqs+1) * sizeof (int *) );
  94.         for(i=0;i<=nseqs;i++)
  95.            sets[i] = (int *)ckalloc( (nseqs+1) * sizeof (int) );
  96.  
  97.         create_sets(0,nseqs);
  98.         fprintf(stdout,"There are %d groups\n",(pint)nsets);
  99.  
  100. /* clear the memory used for the phylogenetic tree */
  101.  
  102.         if (nseqs > 3)
  103.              clear_tree(NULL);
  104.  
  105. /* start the multiple alignments.........  */
  106.  
  107.         fprintf(stdout,"Aligning...\n");
  108.  
  109. /* first pass, align closely related sequences first.... */
  110.  
  111.         ix = 0;
  112.         aligned = (int *)ckalloc( (nseqs+1) * sizeof (int) );
  113.         for (i=0;i<=nseqs;i++) aligned[i] = 0;
  114.  
  115.         for(set=1;set<=nsets;++set)
  116.          {
  117.           entries=0;
  118.           for (i=1;i<=nseqs;i++)
  119.             {
  120.                if ((sets[set][i] != 0) && (maxid[i] > divergence_cutoff))
  121.                  {
  122.                     entries++;
  123.                     if  (aligned[i] == 0)
  124.                        {
  125.                           if (output_order==INPUT)
  126.                             {
  127.                               ++ix;
  128.                               output_index[i] = i;
  129.                             }
  130.                           else output_index[++ix] = i;
  131.                           aligned[i] = 1;
  132.                        }
  133.                  }
  134.             }
  135.  
  136.           score = prfalign(sets[set], aligned);
  137.  
  138. /* negative score means fatal error... exit now!  */
  139.  
  140.           if (score < 0) 
  141.              {
  142.                 return(-1);
  143.              }
  144.           if ((entries > 0) && (score > 0))
  145.              fprintf(stdout,"Group %d: Sequences:%4d      Score:%d\n",
  146.              (pint)set,(pint)entries,(pint)score);
  147.           else
  148.              fprintf(stdout,"Group %d:                     Delayed\n",
  149.              (pint)set);
  150.         }
  151.  
  152.         for (i=0;i<=nseqs;i++)
  153.           ckfree((void *)sets[i]);
  154.         ckfree(sets);
  155.      }
  156.    else
  157.      {
  158. /* clear the memory used for the phylogenetic tree */
  159.  
  160.         if (nseqs > 3)
  161.              clear_tree(NULL);
  162.  
  163.         aligned = (int *)ckalloc( (nseqs+1) * sizeof (int) );
  164.         ix = 0;
  165.         for (i=1;i<=istart+1;i++)
  166.          {
  167.            aligned[i] = 1;
  168.            ++ix;
  169.            output_index[i] = i;
  170.          }
  171.         for (i=istart+2;i<=nseqs;i++) aligned[i] = 0;
  172.      }
  173.  
  174. /* second pass - align remaining, more divergent sequences..... */
  175.  
  176. /* if not all sequences were aligned, for each unaligned sequence,
  177.    find it's closest pair amongst the aligned sequences.  */
  178.  
  179.     group = (int *)ckalloc( (nseqs+1) * sizeof (int));
  180.  
  181.     while (ix < nseqs)
  182.       {
  183.         if (ix > 0) 
  184.           {
  185.              for (i=1;i<=nseqs;i++) {
  186.                 if (aligned[i] == 0)
  187.                   {
  188.                      maxid[i] = 0.0;
  189.                      for (j=1;j<=nseqs;j++) 
  190.                         if ((maxid[i] < tmat[i][j]) && (aligned[j] != 0))
  191.                             maxid[i] = tmat[i][j];
  192.                   }
  193.               }
  194.           }
  195.  
  196. /* find the most closely related sequence to those already aligned */
  197.  
  198.          max = 0;
  199.          for (i=1;i<=nseqs;i++)
  200.            {
  201.              if ((aligned[i] == 0) && (maxid[i] > max))
  202.                {
  203.                   max = maxid[i];
  204.                   iseq = i;
  205.                }
  206.            }
  207.  
  208. /* align this sequence to the existing alignment */
  209.  
  210.          entries = 0;
  211.          for (j=1;j<=nseqs;j++)
  212.            if (aligned[j] != 0)
  213.               {
  214.                  group[j] = 1;
  215.                  entries++;
  216.               }
  217.            else if (iseq==j)
  218.               {
  219.                  group[j] = 2;
  220.                  entries++;
  221.               }
  222.          aligned[iseq] = 1;
  223.          score = prfalign(group, aligned);
  224.          fprintf(stdout,"Sequence:%d     Score:%d\n",(pint)iseq,(pint)score);
  225.          if (output_order == INPUT)
  226.           {
  227.             ++ix;
  228.             output_index[iseq] = iseq;
  229.           }
  230.          else
  231.             output_index[++ix] = iseq;
  232.       }
  233.  
  234.    ckfree((void *)group);
  235.    ckfree((void *)seq_weight);
  236.    ckfree((void *)aligned);
  237.    ckfree((void *)maxid);
  238.  
  239.    aln_score();
  240.  
  241. /* make the rest (output stuff) into routine clustal_out in file amenu.c */
  242.  
  243.    return(nseqs);
  244.  
  245. }
  246.  
  247. int palign1(void)  /* a profile alignment */
  248. {
  249.    int         i,j,k,temp,status;
  250.    int         entries,score,ptr;
  251.    int         *aligned, *group;
  252.    double       dscore;
  253.  
  254.    fprintf(stdout,"\nStart of Initial Alignment\n");
  255.  
  256. /* calculate sequence weights according to branch lengths of the tree -
  257.    weights in global variable seq_weight normalised to sum to 1000 */
  258.  
  259.    seq_weight = (int *) ckalloc( (nseqs) * sizeof(int) );
  260.  
  261.    temp = 1000/nseqs;
  262.    for (i=0;i<nseqs;i++) seq_weight[i] = temp;
  263.  
  264.    distance_tree = FALSE;
  265.  
  266. /* do the initial alignment.........  */
  267.  
  268.    group = (int *)ckalloc( (nseqs+1) * sizeof (int));
  269.  
  270.    for(i=1; i<=profile1_nseqs; ++i)
  271.          group[i] = 1;
  272.    for(i=profile1_nseqs+1; i<=nseqs; ++i)
  273.          group[i] = 2;
  274.    entries = nseqs;
  275.  
  276.    aligned = (int *)ckalloc( (nseqs+1) * sizeof (int) );
  277.    for (i=1;i<=nseqs;i++) aligned[i] = 1;
  278.  
  279.    score = prfalign(group, aligned);
  280.    fprintf(stdout,"Sequences:%d      Score:%d\n",(pint)entries,(pint)score);
  281.    ckfree((void *)group);
  282.    ckfree((void *)seq_weight);
  283.    ckfree((void *)aligned);
  284.  
  285.    for (i=1;i<=nseqs;i++) {
  286.      for (j=i+1;j<=nseqs;j++) {
  287.        dscore = countid(i,j);
  288.        tmat[i][j] = (100.0 - dscore)/100.0;
  289.        tmat[j][i] = tmat[i][j];
  290.      }
  291.    }
  292.  
  293.    return(nseqs);
  294. }
  295.  
  296. double countid(int s1, int s2)
  297. {
  298.    char c1,c2;
  299.    int i;
  300.    int count,total;
  301.    double score;
  302.  
  303.    count = total = 0;
  304.    for (i=1;i<=seqlen_array[s1];i++) {
  305.      c1 = seq_array[s1][i];
  306.      c2 = seq_array[s2][i];
  307.      if ((c1>=0) && (c1<=max_aa)) {
  308.        total++;
  309.        if (c1 == c2) count++;
  310.      }
  311.  
  312.    }
  313.  
  314.    score = 100.0 * (double)count / (double)total;
  315.    return(score);
  316.  
  317. }
  318. int palign2(void)  /* a profile alignment */
  319. {
  320.    int         i,j,k,status;
  321.    int         entries,score,ptr;
  322.    int         *aligned, *group;
  323.  
  324.    fprintf(stdout,"\nStart of Multiple Alignment\n");
  325.  
  326. /* get the phylogenetic tree from *.ph */
  327.  
  328.    if (nseqs > 3)
  329.      {
  330.         status = read_tree(phylip_phy_tree_name, 0, nseqs);
  331.         if (status == 0) return(0);
  332.      }
  333.  
  334. /* calculate sequence weights according to branch lengths of the tree -
  335.    weights in global variable seq_weight normalised to sum to 100 */
  336.  
  337.    seq_weight = (int *) ckalloc( (nseqs) * sizeof(int) );
  338.  
  339.    calc_seq_weights(0, nseqs);
  340.  
  341. /* recalculate tmat matrix as percent similarity matrix */
  342.  
  343.    status = calc_similarities(nseqs);
  344.    if (status == 0) return(0);
  345.  
  346. /* clear the memory for the phylogenetic tree */
  347.  
  348.    if (nseqs > 3)
  349.      {
  350.         clear_tree(NULL);
  351.      }
  352.  
  353. /* do the alignment.........  */
  354.  
  355.    fprintf(stdout,"Aligning...\n");
  356.  
  357.    group = (int *)ckalloc( (nseqs+1) * sizeof (int));
  358.  
  359.    for(i=1; i<=profile1_nseqs; ++i)
  360.          group[i] = 1;
  361.    for(i=profile1_nseqs+1; i<=nseqs; ++i)
  362.          group[i] = 2;
  363.    entries = nseqs;
  364.  
  365.    aligned = (int *)ckalloc( (nseqs+1) * sizeof (int) );
  366.    for (i=1;i<=nseqs;i++) aligned[i] = 1;
  367.  
  368.    score = prfalign(group, aligned);
  369.    fprintf(stdout,"Sequences:%d      Score:%d\n",(pint)entries,(pint)score);
  370.    ckfree((void *)group);
  371.    ckfree((void *)seq_weight);
  372.    ckfree((void *)aligned);
  373.  
  374. /* DES   output_index = (int *)ckalloc( (nseqs+1) * sizeof (int)); */
  375.    for (i=1;i<=nseqs;i++) output_index[i] = i;
  376.  
  377.    return(nseqs);
  378. }
  379.  
  380.